Two-way coupling of FENE dumbbells with a turbulent shear flow 
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We present numerical studies for finitely extensible nonlinear elastic (FENE) dumbbells which are 
dispersed in a turbulent plane shear flow at moderate Reynolds number. The polymer ensemble is 
described on the mesoscopic level by a set of stochastic ordinary differential equations with Brow- 
nian noise. The dynamics of the Newtonian solvent is determined by the Navier-Stokes equations. 
Momentum transfer of the dumbbells with the solvent is implemented by an additional volume 
forcing term in the Navier-Stokes equations, such that both components of the resulting viscoelas- 
tic fluid are connected by a two-way coupling. The dynamics of the dumbbells is given then by 
Newton's second law of motion including small inertia effects. We investigate the dynamics of the 
' flow for different degrees of dumbbell elasticity and inertia, as given by Weissenberg and Stokes 

| numbers, respectively. For the parameters accessible in our study, the magnitude of the feedback 

of the polymers on the macroscopic properties of turbulence remains small as quantified by the 
global energy budget and the Reynolds stresses. A reduction of the turbulent drag by up to 20% is 
observed for the larger particle inertia. The angular statistics of the dumbbells shows an increasing 
alignment with the mean flow direction for both, increasing elasticity and inertia. This goes in line 
^p >' with a growing asymmetry of the probability density function of the transverse derivative of the 

streamwise turbulent velocity component. We find that dumbbells get stretched preferentially in 
, regions where vortex stretching or bi-axial strain dominate the local dynamics and topology of the 

y— j— | velocity gradient tensor. 

GC 

O ' PACS numbers: 47.27.ek, 83.10.Mj, 83.80.Rs 

43 I. INTRODUCTION 

When a few parts per million in weight of long-chained polymers are added to a turbulent fluid its properties 
change drastically and a significant reduction of turbulent drag is observed. [l| Although the phenomenon is known 
from pipe flow experiments for almost 60 years, 0, Q a complete understanding is still lacking. One reason for 
this circumstance is that the physical processes in a turbulent and dilute polymer solution cover several orders of 
! magnitude in space and time; in other words, we are faced with a real multiscale problem. 0, HJ In case of fully 
developed turbulence, the integral scale L, which measures the extension of largest vortex structures in the flow, 
exceeds the viscous Kolmogorov scale r/K, which stands for the extension of the smallest turbulent eddies, by a factor 
of at least 1000. However, long-chained polymers barely exceed the viscous flow scale even in an almost stretched 
state. Their equilibrium extension as given by the Flory radius i?o is usually by a factor of 100 smaller than ?7k- J3] 
In terms of time scales the situation differs slightly. The viscous Kolmogorov time t v can become smaller than the 
slowest relaxation time r of the macromolecules. Although macroscopic closures can rationalize some issues of drag 
• j-h ! reduction Q, the challenging question remains of how the individual dynamics of numerous polymer chains, which is 
present on sub-Kolmogorov and Kolmogorov scales, adds up to a macroscopic effect at scales r < L as being observed 
5^ | in several experiments. [1, d, [l(| 

The description of dilute polymer solutions relies for most studies on one of the following two models: on one 
side, macroscopic continuum models such as Oldroyd-B or FENE-P models (Til H3 , HH, 0, [l5| include the polymer 
dynamics as an additional additive macroscopic stress field. Only the largest scales I > tjk of the viscoelastic fluid 
are described in its full complexity. Numerical problems arise in connection with the pure hyperbolic character of 
the equation of motion for the polymer stress field, such as the conservation of its positivity (see e.g. Ref. [l6| 
for a detailed discussion). In addition, the coarse graining to the macroscopic polymer stress can lead to deeper 
conceptional difficulties, e.g., the failure of energy stability of viscoelastic flows, which is an important building block 
for investigations of stability and upper bounds on the dissipation rate in Newtonian flows. [171 ] Further problems 
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arise for the macroscopic description of non-Newtonian fluids in the limits of very low and high frequencies, where 
they should behave as Newtonian fluids and solids, respectively. [Isl [l9l . l20j 

On the other side, Brownian dynamics models [2lL [22, 0, Hj, l25l| describe the polymer chain on a mesoscopic level 
as overdamped coupled oscillators arranged in bead-spring chains. The models include complex conformations of the 
macromolecules and screening effects due to the solvent such as hydrodynamic interaction. [26[ The simplest of such 
mesoscopic models for a polymer chain is a dumbbell where two beads are connected by a spring. The dynamics in 
these models is on scales £ < t]k- This means that the surrounding fluid is spatially smooth and either a steady [22j|, 
a start-up shear flow (23j, or a white-in-time random flow. [27| In a recent work by Davoudi and Schumacher [28], 
numerical studies at the interface of both descriptions were conducted by combining Brownian dynamics simulations 
(BDS) with direct numerical simulations of a turbulent Navier-Stokes shear flow. The simplest mesoscopic model 
with a linear spring force - the Hookcan dumbbell model - was taken there in order to study the stretching of the 
dumbbell as a function of the outer shear rate and the elastic properties of the springs. However, a feedback of the 
polymers on the shear flow was not included in their study. 

In the following, we want to extend these investigations into two directions. Firstly, we will model the macro- 
molecules more realistically as finitely extensible nonlinear elastic (FENE) dumbbells. Secondly, their feedback on 
the shear flow is included via a two-way coupling. The effect of the FENE dumbbells on the statistical fluctuations of 
the velocity and the velocity gradients will be studied. In addition, conformational properties of the dumbbells, such 
as their extension and angular distribution with respect to the mean flow component, will be addressed. The polymer 
feedback results in an additional forcing that has to be added to the right hand side of the Navier-Stokes equations 
for the advecting Newtonian solvent similar to the case of two-phase flows with dispersed particles [H, [30, HH or 
bubbles. [H| We will keep the full dynamic equation of motion for the dumbbells, containing accelerations due to 
elastic, friction and stochastic forces, and cannot neglect inertia. This step is necessary in order to describe the 
momentum transfer of the dumbbells to the solvent as discussed in Ref. [331 ] . 

In contrast to the conventional BDS that neglect inertia effects from beginning, we will be left here with three 
physical parameters: the Stokes number St for the particle inertia, the Weissenberg number Wi for the elastic 
properties of the dumbbells, and the Reynolds number Re of the flow, respectively The Reynolds number is defined 
as 

Rc = — , (1) 

V 

with the characteristic (large-scale) velocity U, the characteristic length L (both are specified later in the text), and 
the kinematic viscosity of the Newtonian solvent v. The Weissenberg number Wi compares the characteristic dumbbell 
relaxation time r from a stretched to a coiled state with the characteristic time scale of the advecting flow, L/U, and 
is given by 

Wi = ^. (2) 

The Stokes number St relates the particle response time to changes in the surrounding velocity, r st , with the charac- 
teristic flow time scale. It follows to 

St = ^. (3) 

The physics of dispersed FENE dumbbells in a turbulent shear flow is thus described by three dimensionless numbers. 
For a fixed Reynolds numbers Re, we can basically distinguish the following four limiting cases: (i) Wi 3> 1, St ^> 1; 
(ii) Wi < 1, St > 1; (iii) Wi < 1, St < 1; (iv) Wi > 1, St < 1. Case (i) would stand for very heavy particles (or 
dumbbells) which are stretched almost to their contour length. They will behave as dispersed rods. In case (ii), the 
dumbbells would act as heavy spherical particles since they remain coiled in practical terms. The cases of interest for 
dilute polymer solutions are (iii) and (iv), respectively. Inertia effects are then very small, [13] and the Weissenberg 
number can vary from very small to large values implying an increasingly slower relaxation of the macromolecules 
from a stretched non-equilibrium to a coiled equilibrium state in comparison to the characteristic flow variation time 
scale. As we will discuss in the next section, the numerical treatment becomes challenging, on one hand due to the 
finite extensibility, on the other hand due to the small Stokes numbers we are aiming at. The Stokes time r s t sets 
a small but finite time scale then, which can cause stiffness problems for an explicit integration algorithm. Despite 
these efforts, our values for the Stokes number will still exceed the realistic magnitudes for polymer chains in solution 
by orders of magnitude. Nevertheless, we think it is interesting and to some degree necessary to study the dumbbell 
dynamics under these circumstances and to provide a systematic study of how a shear flow will be affected by the 
presence of dispersed bead-spring chains with variable degree of inertia. This will shed some light on possible reasons 
for drag reduction in our model. 
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The outline of the manuscript is as follows. In the next section the equations of motion, the two-way coupling 
and the numerical scheme are presented. Afterwards, we discuss the results for the macroscopic energy balance as 
well as for the Reynolds stresses. This is followed by studies of small-scale properties such as the statistics of the 
extension and orientation of the dumbbells and of their impact on the fluctuations of velocity gradients. We conclude 
with a discussion of our results and will give a brief outlook to extensions of the present work toward more realistic 
parameter settings. 



II. MODEL AND EQUATIONS 
A. The Newtonian solvent 



The Navier-Stokes equations that describe the dynamics of the three-dimensional incompressible Newtonian fluid 
are solved by a pseudo-spectral method using a second-order predictor-corrector scheme for advancement in time.[28| 
The equations of motion are 

— + («-V)u = -Vp + *A7 2 M + / + / P , (4) 
V-u = 0, (5) 

where u is the (total) velocity field, p the kinematic pressure field, / the volume forcing which sustains the turbulence, 
and f p the feedback of the dumbbells (see section II C). The shear flow is modeled in a volume with free-slip boundary 
conditions in the shear direction y and periodic boundaries in the streamwise and spanwise directions x and z. The 
free-slip boundary conditions at y = 0, L y are given by 

«y = °. -dy- = -8y- =0 - < 6 > 

Here, the total velocity field follows by a Reynolds (de)composition as a linear mean part with the constant shear 
rate S and a turbulent fluctuating part 

u = (u) + u' = Sye x + u' . (7) 

The notation (•) stands for the ensemble average, which will be a combination of volume and time averages for most 
cases. The aspect ratio is L x :L y :L z = An: 2: 2ir. The characteristic length is the halfwidth of the slab, L = L y /2. 
Velocities are measured in units of the laminar flow profile U(y) = —\f2cos(ny/2)e x . We will take U x {L y /A) as the 
characteristic velocity U (see also ((T|), ([2J, and ((3|)). The applied volume forcing sustains this laminar flow profile 
and follows from ([4} consequently to f(y) = — \/2ir 2 / (Av) cos(ny /2)e x . Forcing amplitude and profile will remain 
unchanged throughout this study. At sufficiently large Reynolds numbers this linearly stable laminar shear flow 
becomes turbulent when a finite perturbation is applied. 35J The volume forcing / is then a permanent source of 
kinetic energy injection into the shear flow which sustains turbulence in a statistically stationary state. Although 
the steady forcing is of cosine shape, the resulting mean turbulent flow profile will be linear except for small layers 
in the vicinity of both free-slip planes, where the boundary conditions have to be satisfied. Our mean profiles follow 
to (u x (y)) ~ S(y — 1) for y G [0,2] with S = 0.035 — 0.04 for Re = 800. This range of ^-values remained nearly 
unchanged for all parameter sets. In addition, (u' y ) — (u' z ) = 0. The shear flow can be considered therefore as being 
nearly homogeneous. 

The simulation program is run with two spectral resolutions. For Re = 400, a grid with 64 x 32 x 32 mesh 
points was taken. For Re = 800, we took a grid with 128 x 32 x 64 points. The spectral resolution as given by 
the product /c max ?7K = V8ttN x / \2>L x )rj]<^ was 1.5 for the first case and 2.3 for the second. Here, tjk is the viscous 
Kolmogorov scale and defined as t]k = v 3 / 4 / {e') 1 ^ 4 with the mean turbulent energy dissipation rate (e'), where 
e'(x,t) = (v /2)(du' i /dxj + du'^/dxi) 2 for i,j = x,y, z. Clearly, the spectral resolutions are not very large, but they 
give us the opportunity to perform parametric studies in the three-dimensional space which is spanned by Re, Wi, 
and St. Most of our following studies will be conducted for the better resolved case of Re = 800. 



B. The FENE dumbbells 



The smallest building block for the mesoscopic description of the polymer stretching can be accomplished by 
considering dumbbells where two beads (that stand for several hundreds of monomers) are connected by a spring. 
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The entropic elastic force follows the Warner force law [ll| and depends on the separation vector R(t) — X2(i) —Xi(t) 
that is spanned between both beads at positions x 2 (t) and Xxft), respectively. The force law is given by 

Fci{R) = i- H */Lr (8) 

where Lq is the contour length of the dumbbells which cannot be exceeded. The spring constant is denoted by H . When 
taking into account the elastic entropic force, hydrodynamic Stokes drag, and thermal noise, the second Newtonian 
law for a FENE dumbbell written in relative coordinates R{t) and center-of-mass coordinates r(t) = (xi(t) + X2(t))/2 
reads [13, H3 

r = v, (9) 
Y v = -v + -(u 1+ u 2 ) + J^Z r , (10) 
R = V, (11) 



m b ■ xr A 2HR Uk B T „ 

T V.-V + A»- al _ R2/Lg) +v /— fa, (12) 

where Am = u(x2,t) — u(xi,t) is the relative fluid velocity at the bead centers. The last terms in the velocity 
equations, containing and £r, stand for vectors of thermal Gaussian noise with the properties 

(&(*)) = 0, (13) 
= 5 l3 5(t-t') (14) 

for i, j = x,y, z. The three components of each vectorial noise term are statistically independent stochastic processes. 
Furthermore, the vectorial noise with respect to the center-of-mass velocity is statistically independent to that for the 
relative velocity dynamics. The noise prevents the extension of a dumbbell to shrink below its equilibrium length 



with fee being the Boltzmann constant, T the temperature. Equation (I15| follows from the equipartition theorem. 
The contour length Lq = 10R is used throughout this study and Rq — rj^. The relaxation time of the dumbbells is 
given by [III 

where 

( = 6np[va (17) 

is the Stokes drag coefficient of a spherical bead with radius a. The fluid mass density is pf. Due to the current 
resolution contraints the dumbbells will experience both the smooth and partly rough scales of the advecting flow. 
Consequently, the velocity difference Am is kept in the equation and not approximated by the linearization Am m 
(R ■ V)m as it is done in BDS where L <C t\k ■ For spatially smooth flows both expressions give the same results. 

The equations © through (TT!?|) introduce the other two dimensionless parameters beside the Reynolds number Re, 
the Weissenberg number Wi and the Stokes number St, respectively (see definitions §2$ and ©). The Stokes time r s t 
is the response time of an inertial particle which is required to speed up to the velocity of its local surrounding. A 
zero Stokes time implies a behavior as a passive Lagrangian tracer. For beads, this time follows to r st = mb/C with 
C as given above and consequently 

T st = . (18) 

The density contrast p p /pt is to very good approximation unity 36], i.e. polymers are considered as neutrally buoyant. 
In Ref. [28j |. we have compared the polymer relaxation time to the microscopic stretching time scale. This is given 
by the inverse of the maximum Lyapunov exponent and is comparable to the microscopic time scale of the flow, the 
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Wi = 


= 3 




Wi„ 


= 0.8 




Wi„ 


= 0.6 




Wi = 


= 20 
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= 5.1 




Wi„ 


= 4.3 




Wi = 


= 100 




Wi„ 


= 25.7 




Wi„ 


= 21.5 




St = 


5.0 x 10" 


-4 


St^ - 


= 1.3 x 10" 


-4 


St^ - 


= 1.1 x 10" 


-4 


St = 


5.0 x 10" 


-3 


St-q - 


= 1.3 x 10" 


-3 


St^ - 


= 1.1 x 10" 


3 


St = 


5.0 x 10" 


-2 


Stjj - 


= 1.3 x 10" 


-2 


St^ - 


= 1.1 x 10" 


-2 


St = 


5.0 x 10" 


-1 


St^ - 


= 1.3 x 10" 


-1 


St^ - 


= 1.1 x 10" 


-1 



TABLE I: The Weissenberg and Stokes numbers rescaled by the Kolmogorov time t v of the flow. Wi,, = t/t v and St,, = r s t/r, 
Note that t v is based on the pure Newtonian case. Only minor changes arise when polymers are added to the solvent. 



Kolmogorov time t v = y/i//(e). Table 1 gives an overview of the values of St and Wi that have been used and of how 
they translate into St,, and Wi,,, respectively. We see that the Stokes numbers get as low as 10~ 4 when measured in 
viscous units, which is still orders of magnitude above the realistic estimates for dilute polymer solutions which are 
about three to four order of magnitude below our minimal value. 

In most cases, an ensemble of 6.3 x 10 FENE dumbbells, i.e. 1.2 x 10 5 beads, is advanced by a weak second-order 
predictor-corrector scheme simultaneously with the flow equations. (2~0| The finite extensibility and the small Stokes 
numbers require a semi-implicit time-stepping for some variables. In order to avoid a total length larger than Lo, 
we proceed in line with Ref. [IjJ and solve a cubic equation for R = \R\ in the corrector step. Initially, the center 
of mass of the dumbbells is seeded randomly in space with a uniform distribution and an initial extension of Rq. 
All Lagrangian interpolations were done with a trilinear scheme. Details on the numerical procedure are outlined in 
appendix A. 

In order to build a bridge to macroscopic simulations we provide an estimate for the contribution of the dumbbell 
ensemble to the zero-shear viscosity. Following Ref. [2l[ it is defined as 

V P = Pp v p = n p k B Tr , (19) 

with the number density of dumbbells n p . When applying (fT5|) as well as definitions (fl"6]) and (JTTj) , and using pf/p p = 1 
one gets 

3 

v p = --KUpRlua (20) 

with the solvent viscosity v. The bead radius a is substituted by the Stokes time r st . Recalling the definitions for the 
Kolmogorov length ??k = v 3 ^ 4 / (e') 1 / 4 and for the Kolmogorov time t v = y/v/{e'), one ends with the relative viscosity 

S = ~v = ^"p^V^- ( 21 ) 

For the present simulations, one dumbbell is seeded per grid cell and therefore n p ~ Additionally, Rq ~ t)k- 

Following table 1 for the runs at Re = 800, one gets ratios of s between between 0.1 for the smallest Stokes number 
going up to 3 for the largest one. The latter value is rather large for polymer solutions. Values below unity are usually 
taken, such as in DNS with the Oldroyd-B model, [l 4] Equation ([2"Tj) is in this spirit consistent with the discussion 
in the introductory part. Only the lower Stokes numbers result to values of s as taken for macroscopic DNS for 
viscoelastic shear flows. 



C. Two-way coupling 



The back-reaction of the dumbbells on the fluid consists of contributions from the Stokes friction and the stochastic 
noise term. In accordance with Newton's third law, the force contribution from each of the two beads at positions Xi 
(i = 1, 2) follows to 

Ft = -F q {st) - F[ n) = (( Xl - u( Xi )) - y2W;& ■ (22) 
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The force density generated by all FENE dumbbells results to 

N p 2 

/*/p = EE^-*f), (23) 
j=l i=l 

where N p is the number of dumbbells. The volume integral of (f2"5|) gives a force since the delta function carries the 

dimension of an inverse volume due to J S(x — x^)d 3 x = 1. Consequently, the dimensionless form of the forcing 
reads 

6{x-x^), (24) 

where the bead volume follows to Vb = 47ra 3 /3 = (47r/3)(9!/r s t/2) 3 / 2 . The notation 5 is for the dimensionless delta 
function. We have used again pi/p p « 1. The force density has to be evaluated at space points that are between the 
mesh vertices. Again the trilinear interpolation has to be used to evaluate the contributions of the point force to the 
eight next neighboring mesh vertices. 

III. LARGE-SCALE PROPERTIES 
A. Energy balance 

The first analysis step is the study of the effects of the two-way coupling on the macroscopic properties of turbulence. 
Given the boundary conditions for our problem, eq. (U) results in the following balance for the total kinetic energy 

E ( t ) = wJv M 2 d3x with V = L * L y L z, 

dE 

— = -^{(dujdx^v + (u ■ f) v + (u ■ f p ) v , 

= -e(t) + e- m (t) - E p (t) (25) 

where (-)y = y J -d 3 x is the short notation for the volume average. In case of statistical stationarity, one gets 
d(E) t /dt = and thus 

(e in ) = (e) + (e p ) . (26) 

Figure [T] shows the three mean rates as a function of the Stokes number for two Weissenberg numbers Wi = 20, 100. 
The mean energy dissipation rate (e) and the mean energy injection rate (£m) are of the same order of magnitude for 
all cases. They remain nearly unchanged with respect to Weissenberg number, which indicates that the effect of the 
dumbbell ensemble on the macroscopic flow properties is small. Nevertheless, one observes a slight increase of the 
mean energy injection rate (£i n ) with respect to St going in line with a decrease of (e) (see upper and mid panel of 
Fig. [I}. Recall that the energy injection rate will be maximal for the laminar case, i.e. for u \\ f. The trend of the 
data indicates that the streamwise flow component relaminarizes slightly with growing inertia. The lower panel of the 
same figure shows the findings for the dissipation due to polymer stretching (e p ). As an additional energy dissipation 
mechanism, it consumes injected energy which goes into the elastic energy budget of the dumbbell ensemble. The 
rate (e p ) grows in magnitude with respect to both parameters, the Stokes and Weissenberg number. For Wi = 3, the 
dumbbells are not significantly extended and no clear trend of (e p ) with St could be observed. The dissipation rate 
(e p ) is significantly smaller in comparison to the runs with larger Wi. 

In order to estimate the maximum feedback of the dumbbells on the flow, we performed an "academic experiment" 
for our system by tethering one of the two beads of a dumbbell at a fixed position. The dumbbells get then stretched 
more efficiently and undergo strong conformational fluctuations. Figure 2 illustrates their dramatic effect on the total 
kinetic energy. We compare the freely draining case with the tethered one and observe a significant decrease of the 
kinetic energy. An inspection of the flow structures indicates that the turbulent fluctuations are supressed almost 
completely. The flow becomes basically laminar. The magnitude of the feedback for freely draining dumbbells will 
always remain significantly below this artifical limit with tethered dumbbells. 
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B. Reynolds stresses 

Figure [3] shows the four non-vanishing components of the Reynolds stress tensor (u' i u'j)/(2k) where k — ((u' i ) 2 }/2 
is the turbulent kinetic energy (TKE). The moments are averages over the whole simulation volume for a sequence of 
about 100 statistically independent snapshots of the time evolution of the shear flow. The results can be summarized 
to the following trends. For the two smallest Stokes numbers, no dependence on the Weissenberg number is observed. 
For St = 0.05 and 0.5, the mean streamwise fluctuations are enhanced while the remaining components of the Reynolds 
stress tensor decrease as a function of Wi. This finding is in agreement with observations in a Kolmogorov flow by 
Boffetta et al. [37| 

Similar to the friction factor for a turbulent pipe [38j ]. we can define a friction factor for the present flow where 
the applied pressure gradient term has to be substituted by an amplitude of the static volume forcing profile / that 
sustains the laminar cosine flow profile. Consequently, 

IF L 

Cf = {u x (y = l y ))- ■ (2?) 

Since f(y) — — \f2~7r 2 / (Av) cos(ny /2)e x , we take F — f x (y — L y ) — \f2-K 2 j{Av). A similar definition was suggested for 
a Kolmogorov flow which is also driven by a volume forcing. [37] Drag reduction by dispersed dumbbells would go in 
line with a decrease of the dimensionless measure c/ below the Newtonian value For the smallest Stokes number, 
the ratio goes to about unity. The slight overshoot is attributed to the strong variations of the streamwise velocity 
at the free-slip planes. Figure 0] indicates a reduction by 20% — 25% at St = 0.05,0.5 and for the larger Weissenberg 
numbers. The series with Wi = 3 gave c/ ~ . 

An important structural ingredient of shear flows are the asymmetric fluctuations of the three diagonal elements 
of the Reynolds stress tensor. The streamwise fluctuations ((u' x ) 2 ) are spatially arranged in streamwise streaks which 
interact with streamwise vortices in a so-called regeneration cycle of coherent structures. This cycle is sustained 
by the non- normal amplification mechanism. [39l l4fj| The impact of long-chained polymers on the extension of the 
streamwise streaks has been demonstrated in experiments 10] and numerical simulations. [4l|, |42j While streamwise 
fluctuations were found to increase, the fluctuations in shear and spanwise directions decreased. This is in line with 
our observations as discussed above. In Fig. [5J we show isolevels of the streamwise turbulent fluctuations for opposite 
sign at Wi = 3, 20, 100. Although not very pronounced, a slight increase in the connectivity and extension of the 
streamwise streaks can be observed with increasing Weissenberg number. 

As we can see, the statistics of macroscopic turbulent properties is affected only slightly by the dispersed FENE- 
dumbbells. Their impact increases with Weissenberg number as well as with Stokes number. In order to rule out 
that particle inertia dominates the discussed trends of our studies, we considered the case of dispersed beads in the 
same flow at the same Stokes numbers. This is achieved by switching off the elastic spring force, i.e. F c \ = 0. The 
Stokes friction force remained as the only force. The quantity f p models then the feedback of the particles on the 
flow. We added the statistical means of injection and dissipation rates as a function of the Stokes number for this 
case to Fig. 1. While the mean injection and mean dissipation rates are of the same magnitude, the dissipation due 
to particle feedback is orders of magnitude smaller in comparison to the polymer feedback, except for the largest St. 
In addition, we found no clear trends for the Reynolds stress components as a function of St. 



IV. SMALL-SCALE PROPERTIES 



A. Extensional and angular statistics of dumbbells 

The finite extensibility of the dumbbells will affect the shape of the probability density function (PDF) of R, which 
is supported on scales smaller than Lq only. Figure [5] reports our findings for p(R) for different Weissenberg and 
Stokes numbers. For the lowest Weissenberg number, Wi = 3, the majority of the dumbbells remains at the extension 
of about the Kolmogorov length t]k- This picture changes for larger values of Wi. At Wi = 100, the majority of the 
ensemble is stretched to almost Lq, which manifests in the sharp maximum at R SLq- Qualitatively, the ch ang e of 
the shapes of the PDFs with increasing Wi agrees well with experimental findings [43| and analytical studies [U |44[ 
for the coil-stretch transition in random flows. The trends with the Stokes number remain small in all cases. However, 
the data show that growing particle inertia suppresses the stretching to very extended molecules since the response 
time of the molecules to the variation of the structures increases (see e.g. mid panel of Fig. |fj]). 

As we have seen in the last section, the fluctuations of the turbulent velocity field in the shear flow vary strongly from 
one space direction to another (see e.g. Fig. [3]). The major contribution is contained in the streamwise component 
{(u' x ) 2 } parallel to the direction of the mean turbulent flow. This suggests an investigation of the angular statistics of 
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the polymers since their stretching can be expected to become anisotropic as well. The following dumbbell coordinate 
system will be used therefore throughout this text: R x = R cos p cos 9, R y = R sin p cos 9, and R z — RsinO, where 
R is the distance between both beads. The notation differs from conventional spherical coordinates, but has the 
advantage of giving perfect alignment with the outer mean flow direction for ip = 9 = 0. (p is the azimuthal angle and 
9 the polar angle. While the azimuthal angle always remains in the shear plane that is spanned by the streamwise 
and shear directions, the polar angle 9^0 indicates a dumbbell orientation out of this plane. 

Davoudi and Schumacher [28| discussed the statistics of both angles as a function of the Weissenberg number for 
passively advected Hookean dumbbells. The PDF of the polar angle was found to remain symmetric and to be less 
sensitive with respect to variations of Wi. Our focus will be therefore on the statistics of the azimuthal angle ip 
which can take values between —tt/2 and tt/2. The asymmetry between both quadrants is quantified by the following 
measure for the PDF p(ip): 

A(<p)=p(<p)-p(-p), (28) 

with ip € [0,7r/2]. The measure A(p) is plotted for two Weissenberg numbers in Fig. [71 A pronounced maximum 
of A{p) implies that the dumbbells are preferentially slightly tilted in the direction of shear, away from the mean 
flow direction (see an illustration in Fig. [5]) . We find that with increasing Weissenberg number the asymmetry of the 
angular distribution grows in magnitude. The same trend holds when the Stokes number grows at fixed Weissenberg 
number. In each case, the graph of A(p) shows an increasingly sharper maximum, which is shifted towards smaller 
p. Fluctuations of the dumbbells in the vicinity of ip — are enhanced while the tails for very large p are depleted. 
Growing inertia amplifies this trend. Once the dumbbells are aligned along the mean flow they remain in this 
orientation for longer periods of their evolution. 



B. Velocity gradient statistics 



Since the polymer dynamics takes place at the smallest scales of the turbulent flow, we study the impact of the 
dumbbells on the small-scale statistical properties of the flow in the following. Recent experimental and numerical 
studies in simple Newtonian shear flows indicate that in particular the statistics of the transverse derivative of the 
streamwise turbulent velocity component du' x /dy is a sensitive measure for detecting deviations from local isotropy 
in homogeneous or nearly homogeneous shear flows. [IH H|| |47} In a shear flow with a mean shear rate S > 0, one 
expects a positive value for derivative skewness and other higher odd order moments which are defined as 

^^> = J|» 

The derivative moments would be exactly zero in a perfectly isotropic flow. Their non-zero magnitudes indicate that 
velocity gradient fluctuations of the streamwise component along the direction of the outer shear gradient are more 
probable than the ones in the opposite direction. It can be expected that the asymmetry in the angular distribution, 
which we discussed above, will have an impact on the statistics of exactly these gradient fluctuations. Figure [5] 
reports our findings for the PDF of the transverse derivative, which has been normalized by its root mean square 
value for all cases. We observe in both figures a depiction of the left hand tail, which stands exactly for the velocity 
gradient fluctuations opposite to the direction of the mean shear. The results suggest that the preferential orientation 
fluctuations of the dumbbells at azimuthal angles p > go in line with a depletion of the negative tail of the PDF 
of the transverse derivative. As sketched in Fig. [8j negative transverse gradients would be amplified by prefential 
orientations with p < which correspond to the dumbbell colored in gray. The findings are consistent with our 
observations on the ^-statistics . They can also be rationalized (but not explained) when considering the equation for 
the Brownian dynamics of the FENE dumbbell [2l[ 

In the plane shear flow geometry the component R x along the mean flow direction is of particular interest. Since we 
are interested in stretched dumbbells with R x > i? an d in Wi > 1 we neglect contributions from the spring force and 
the noise for a moment. With the Reynolds decomposition ([7J we get 

TT=( s + £)* + (g)* + - 
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The important term is the first term on the r.h.s. of (|3Tj) . The other three contributions will behave as noise terms. 
Fluctuating gradients du' x /dy along Se y lead to a more rapid growth of R x (for an angle ip > 0) and a prefered 
alignment with the mean flow. This causes a more rapid decrease of R y and consequently of R x via (|31jl . The dumbell 
can be kicked afterwards again to larger ip values and transfers momentum to the flow which corresponds exactly to a 
local patch of du' x /dy > (see also Fig. ©). Then R y grows and this whole cycle starts anew. Small scale gradients 
with the opposite sign diminish the total shear in the surrounding of the dumbbell and cause a less efficient stretching 
and cycle. Clearly, this picture omits some important features such as the tumbling of the dumbbells. 

The depletion of gradient fluctuations goes in line with experimental observations by Liberzon et al. p8l . |49] | The 
authors found e.g. that the enstrophy production became anisotropic when polymers are added to the fluid. This 
quantity is directly related to transverse gradient components discussed here. 



C. Invariants of the velocity gradient tensor and dumbbell extension 

The efficient stretching of the dumbbells is connected to particular local flow topologies. They are related to the 
three eigenvalues Ai of the velocity gradient tensor or the corresponding three velocity gradient tensor invariants, 
which are denoted as I\, I2, and I3. The eigenvalues of the velocity gradient tensor du^/dxj result as zeros of the 
following third-order characteristic polynomial 50] 



For an incompressible flow [5S 



A 3 -/iA 2 +/ 2 A-/ 3 = 0. (33) 



h = A 1 + A 2 + A 3 = Tr( ) =0, 



dxj 



A1A2 + A2A3 + A3A1 



1 du' du'; 







2 dxi dxi ' 



du' t \ _1 du'i du'j du' k 



h = AiA 2 A 3 = det ^ = - ^ -± -jj-* . (34) 
\ axj J i oxj axk oxi 

The remaining coefficients of (|33j) are therefore I2 and .Z3, which span the I3 — I2 parameter plane. The scatter plots 
for turbulent flows result in a typical skewed teardrop shape. With our definitions given above the following crude 
classification scheme can be given. For I2 > 0,^3 > vortex stretching is present corresponding to Ai = a, A2,3 = 
—a ± ib (first quadrant); for I2 > 0, 13 < vortex compression is present corresponding to Ai = —a, A2.3 = a ± ib 
(second quadrant). The cases I3 < are associated with bi-axial strain for 7 2 < (third quadrant) corresponding to 
Ai = a, A2 = b, A3 = — (a + b) and with uniaxial strain at I2 > (fourth quadrant) corresponding to Ai = a, A2 = 
—6, A3 = —(a — b). Constants a and b are larger than zero in all cases. Figure [TO] relates the extension of the dumbbells 
to the corresponding local velocity gradients in the I3 — I2 plane (and consequently to the existing local flow topology). 
The invariants of the velocity gradient were evaluated in the center of mass of each dumbbell. The typical teardrop 
shape for the turbulence data in the parameter plane is detected. 

Our findings can be summarized as follows. Strongly stretched dumbbells go in line with the largest excursions of 
the gradients in the ^3 — 1% plane. The longest dumbbells are found preferentially in regions where vortex stretching 
or bi-axial strain dominate the local flow topology. The preferential stretching by bi-axial strain was discussed already 
for the passive advection of FENE dumbbells in a minimal flow unit.[25[ It corresponds to the scenario that different 
parts of the dumbbell get pulled by counterstreaming streamwise streaks. The preferential extension close to vortex 
stretching means that the polymers are pulled around streamwise vortices. This point was outlined in Ref. [42| on 
the basis of an analysis of the energetics of viscoelastic turbulence. Here, we find both in a common description 
based on the analysis of the full velocity gradient tensor, i.e. the symmetric strain tensor plus the anti-symmetric 
vorticity tensor. We do also observe that the area of the teardrop shape shrinks with increasing Stokes number. This 
indicates that the small-scale velocity gradients are supressed in magnitude, which goes in line with more limited 
excursions across the ^3 — I2 plane and a rclaminarization of the turbulence. Again, this goes in line with very recent 
experimental observations by Liberzon et aLpfoT] 



V. SUMMARY AND DISCUSSION 



The presented numerical studies aimed at connecting a macroscopic description for the Newtonian turbulent shear 
flow to the mesoscopic description of an ensemble of FENE dumbbells which are advected in such flow. The momentum 
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transfer of the dumbbells with the fluid is implemented by an additional volume forcing in the Navier-Stokes equations. 
In numerical terms, pseudospectral simulations for the solvent are coupled to a system of stochastic nonlinear ordinary 
equations in order to model a viscoelastic fluid. 

For the accessible parameters we found slight modifications of the macroscopic flow structures and mean statistical 
properties only. This was demonstrated for the global energy balance and the mean components of the Reynolds 
stress tensor. We conclude that dumbbell inertia effects are present, but remain sublcading in comparison to the 
elastic properties. For the present viscoelastic flow a drag reduction of up to 20% is achieved. The microscopic 
properties of turbulence were found to be more sensitive with respect to the Weissenberg number. The statistics of 
the azimuthal angle ip is consistent with former findings for elastic Hookean dumbbells. (28j A growing number of 
dumbbells becomes increasingly aligned with the mean flow direction. The feedback of the FENE dumbbells on the 
small-scale properties of turbulence is demonstrated for two gradient measures, the PDF of the transverse derivative of 
the turbulent streamwise velocity component du' x /dy and the diminished scattering of the velocity gradient invariants 
amplitudes in the ^3 — I2 plane with increasing Wi. The asymmetry of the PDF p(du' x /dy) is found to increase with 
increasing Wi. Furthermore, we determined that strongly stretched dumbbells can be found close to vortex stretching 
or biaxial strain topologies of the advecting shear flow. 

The present study should be considered as a first step for such class of hybrid models. One difference to the 
situation in a dilute polymer solution is the relatively large Stokes number that had to be taken. Our dispersed 
dumbbells behave in parts like deformable particles rather than polymer chains. Frequently heavier quasi-particles 
are used for the study of turbulence in particle-ladden flows. (3l| Extensions of our investigations will have to go into 
two directions. Firstly, it is desirable that larger spectral resolutions, like the ones in Ref. [28[, are achieved. This 
will require a fully parallel implementation of the current numerical scheme. Larger computational grids and higher 
Reynolds numbers will give us the opportunity to decrease the ratio i?o/ 7 7K and to increase Lq/Rq to more realistic 
values. Secondly, eq. (|21|) implies the efforts that have to be taken in order to approach the situation in a polymer 
solution. Decreasing values of Rq and St have to be compensated by n p , e.g., a reduction of both - Rq and St r) - 
by an order of magnitude requires an increase of the concentration (or number density) by a power of 5/2. Once 
such operating point is reached, the time scale argument which is thought to be important for the drag reduction 
effect, can also be studied. [1] Finally, a recent work by Vincenzi and co-workers [511 ] provides an interesting ansatz for 
modelling the polymer dynamics. The authors studied a conformation-dependent Stokes drag coefficient that caused 
a significant dynamical slow-down of the coil-stretch transition in steady elongational and random flows. The test of 
these ideas in turbulent shear flows is still to be done. 
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APPENDIX A: SEMI-IMPLICIT INTEGRATION SCHEME FOR DUMBBELLS 

The FENE dumbbells consist of two beads at positions x\(t) and x%(£) which are connected by a nonlinear elastic 
spring. The velocities of the advecting flow at both beads are denoted by U\ and u%, respectively. Note that these 
velocities coincide with x\ and £2, respectively, for St = only. Since the beads are usually found between mesh 
vertices, the values for U\ and u 2 have to be determined by trilinear interpolation from the known velocity vectors 
at the neighboring grid sites. The dynamical equations for the dumbbells are set up in relative and center-of-mass 
coordinates. The relative coordinate (or separation) vector of the dumbbell is given by 

R(t) =x 2 (t) -xi(t). (Al) 

The center-of-mass coordinate vector is given by 

r(t) = -(x 1 (t)+x 2 (t)). (A2) 

The velocities which are assigned with the relative and center-of mass coordinates are denoted as V and v, respectively. 
The Newtonian equations for the dynamics of the FENE dumbbells in dimensionless form, which follow then from 
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([5])-([T2]l with the definitions {3) and are given by 



df 






= v, 


717 


dv 


1 


dt 


St 


dR 






= v, 


~dT 


dV 


1 


~dT 


St 



1 , _ _ , i?o r 
-t> + -(m+ti 2 ) + — ^= ^r 

2 2LyWi 



-V + («2 - fix) - 



ft 



2Wi 



i(l -R 2 L 2 /L 



<WiL 



(A3) 
(A4) 

(A5) 
(A6) 



For the following, we omit the tilde symbol for the dimensionless quantities. The predictor values of the center-of-mass 
vector r and the distance vector R are calculated by an explicit Euler step whereas the corresponding velocities are 
treated by an implicit Euler step, giving 



r l + Atv 
1 



v = 



St + At 
R* = R l + AtV l , 
1 



St v l + -(u\ + u l 2 )At + R ° Aw l 
2 V 2LVWI 



V* = 



St + At 



St V' + K - «i)At - 



R< 



2Wi(l- {R l ) 2 L 2 /Ll 



At- 



Rq 



LVWi 



:AW l 



The corrector step for the center-of-mass and distance vectors is given as 



J+i 



..'+1 



r l + i(t>* +v l )At 



(A7) 
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(All) 
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(A12) 
(A13) 
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- I St V* + (u* 2 - w*)At + St V 1 + (u\ - u[)At 



R l+1 
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R 1 
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At 



Rq 



LVWi 



-.AW 1 



(A14) 



Note that the corrector step for the distance vector is semi-implicit in the velocity in order to avoid stiffness of the 
equation system at small Stokes numbers. The corrector step for the distance velocity V has to be semi-implicit in 
the separation vector R due to the finite extensibility of the dumbbells. [2~fl] When inserting (|A14[) into (|A13[) one gets 



1 + 



(Atf 



8Wi (St + At) (1 - ( J R'+ 1 ) 2 L 2 /Lg 



R 



i+i 



(A15) 



where the abbrevation A contains terms only which are known. By taking the norm of (|A15|) one ends up with a 
cubic polynomial for R l+1 . The formula for the "casus irreducibilis" of three real solutions of the polynomial goes 
back to F. Viete [52j and yields directly the unique solution for R = \R\ between and Lq. From (|A15|1 follows now 



R i+i = R i+i 



A 



(A16) 
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This value is inserted into (|A14[) which completes the corrector step. 
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FIG. 1: Mean dissipation and injection rates as a function of the Stokes and Weissenberg numbers. Upper picture: mean energy 
injection rate due to shear flow forcing (gin)- Mid panel: mean energy dissipation rate (e). Lower panel: mean dissipation rate 
which arises from the coupling to the dumbbell ensemble (e p ). The Reynolds number is Re = 800. The case with F e \ = is for 
the case without spring force and stands for a shear flow with dispersed inertial particles. 
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numbers are shown. The Stokes numbers of the data are indicated in the legend. Data are for Re = 800. 
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FIG. 7: Asymmetry of the probability density function (PDF) of the azimuthal angle ip. It is defined as A((p) = p(f) — p(—ip). 
The upper panel shows the data for Wi = 3 and four different Stokes numbers. The lower panel shows the data for Wi = 100 
and four different Stokes numbers. The analysis is for Re = 800. 
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FIG. 8: Sketch of the orientation of a dumbbell in the turbulent shear flow. The mean turbulent flow profile is indicated. 
The dark-colored dumbbell stands for the preferentially oriented one while the gray-colored orientation is less probable. This 
orientation asymmetry leads to the asymmetry in the angular distribution as given in Fig. [7] 
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FIG. 9: Probability density function (PDF) of the transverse velocity gradient of the streamwise turbulent fluctuations, du' x /dy. 
The Newtonian case is compared with the two larger values of the Weissenberg number at St = 5 x 10~ 4 and 0.5, respectively. 
The data are for Re = 800. 
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FIG. 10: Relation between the extension of the dumbbells and the local velocity gradient at the center of mass of the dumbbells. 
The local flow topology that is related to the velocity gradient is quantified by the second and third invariants I2 and 73 (see 
eqns. I|34|) for the definition). Quadrant I stands for vortex stretching, II for vortex compression, III for bi-axial strain, and IV 
for uniaxial strain, respectively. The gray color coding of the bins for < R/Lo < 0.25, 0.25 < R/Lo < 0.5, 0.5 < R/Lo < 0.75, 
0.75 < R/Lq < 1 is indicated by the legend for each figure. Data are for Re = 800. 



